Preamble

library(doParallel)
## Loading required package: foreach
## Loading required package: iterators
## Loading required package: parallel
library(tictoc)
library(rstan)
## Loading required package: StanHeaders
## Loading required package: ggplot2
## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
library(RColorBrewer)
library(survival)
doParallel::registerDoParallel(cores = 8)
FORCE_RE_RUN = F

# set of functions for the simulations
source('functions.R')
# get posterior samples from the pharmacodynamic model
load('Rout/stan_out_rate.RData')

Plot model output

Figure 1

### Figure to demonstrate what the model looks like
sim_dat_null = simulate_data(N = 10^4, effect = 1, thetas = thetas, xs = 0:21)
sim_dat_5 = simulate_data(N = 10^4, effect = 1.05, thetas = thetas, xs = 0:21)
sim_dat_75 = simulate_data(N = 10^4, effect = 1.075, thetas = thetas, xs = 0:21)
sim_dat_15 = simulate_data(N = 10^4, effect = 1.15, thetas = thetas, xs = 0:21)
sim_dat_25 = simulate_data(N = 10^4, effect = 1.25, thetas = thetas, xs = 0:21)


par(las = 1, mfrow=c(2,2), family='serif', bty='n',cex.axis=1.4,cex.lab=1.4)
set.seed(485)
plot(0:21, sim_dat_null[1,], col = adjustcolor('grey',.3),type='l',ylim = c(40,13),
     ylab='Ct',xlab = 'Days since peak viral load')
points(0:21, sim_dat_null[1,], pch='.')
cols = RColorBrewer::brewer.pal(n = 11, name = 'Set3'); j = 1
for(i in sample(0:1000, 9)){
  lines(0:21, sim_dat_null[i,], col = adjustcolor(cols[j%%11 + 1],.4))
  points(0:21, sim_dat_null[i,], pch='.')
  j = j+1
}
mtext(text = 'A', side = 3, adj = 0, cex=2)

cols = rev(RColorBrewer::brewer.pal(n = 5, name = 'RdYlBu'))
plot(0:21, colMeans(sim_dat_null), type='l', ylab='Ct', col=cols[1],
     xlab = 'Days since peak viral load', ylim = c(40,13),lwd=2)
lines(0:21, colMeans(sim_dat_5),col=cols[2],lwd=2)
lines(0:21, colMeans(sim_dat_75),col=cols[3],lwd=2)
lines(0:21, colMeans(sim_dat_15),col=cols[4],lwd=2)
lines(0:21, colMeans(sim_dat_25),col=cols[5],lwd=2)
legend('topright', legend = c(0,5,7.5,15,25),col=cols, bty='n',
       inset=0.03,lwd=3, title = 'Effect size (%)')
mtext(text = 'B', side = 3, adj = 0, cex=2)


rates_15 = compute_rate(sim_data = sim_dat_15, xs = 0:21)
rates_0 = compute_rate(sim_data = sim_dat_null, xs = 0:21)

plot(density(rates_0), main = '', ylab = '', col = cols[1],lwd=3,
     yaxt='n', xaxt='n', xlab='Clearance half-life (hours)', xlim = c(3,1))
lines(density(rates_15), col=cols[4], lwd=3)
axis(1, at = seq(1,3,by=.5), labels = 24/seq(1,3,by=.5))
legend('topright', legend = c(0,15),col=cols[c(1,4)], bty='n',
       inset=0.03,lwd=3, title = 'Effect size (%)')
mtext(text = 'C', side = 3, adj = 0, cex=2)


TC_0 = compute_clearance_time(sim_dat_null,xs = 0:21)[,1]
TC_15 = compute_clearance_time(sim_dat_15,xs = 0:21)[,1]
plot(density(TC_15[complete.cases(TC_15)], bw = .6),
     main = '', ylab = '', col = cols[4],lwd=3,
     yaxt='n',  xlab = 'Clearance time (days)')
lines(density(TC_0[complete.cases(TC_0)],bw = 0.6), col=cols[1], lwd=3)
legend('topright', legend = c(0,15),col=cols[c(1,4)], bty='n',
       inset=0.03,lwd=3, title = 'Effect size (%)')
mtext(text = 'D', side = 3, adj = 0, cex=2)

Set Parameters

# Parameters for the simulations
Sample_sizes = seq(20,100,by=10)
Effect_sizes = c(1.05,1.075,1.1,1.15,1.2,1.25)
Max_follows = c(7,10,14) # Follow-up durations
Nsim = 2000 # number of trials
my_alpha = 0.05 # significance level to compute type 2 error

Simulation 1

Compare power for a bunch of sample sizes, effect sizes and differing follow-up durations. Takes about 1-2 hours on 8 cores.

out_path = 'Rout/power_summaries.RData'
if(FORCE_RE_RUN | !file.exists(out_path)){
  power_summaries = list(); j = 1
  for(i in 1:length(Max_follows)){
    pp1 = compute_pvalues(Effect_sizes = Effect_sizes,
                          Sample_sizes = Sample_sizes,
                          Follow_up_days = 0:Max_follows[i],
                          Nsim = Nsim,
                          summary_function = compute_clearance_time,
                          thetas=thetas, endpoint = 'timetoevent')
    pp2 = compute_pvalues(Effect_sizes = Effect_sizes,
                          Sample_sizes = Sample_sizes,
                          Follow_up_days = 0:Max_follows[i],
                          Nsim = Nsim,
                          summary_function = compute_rate,
                          thetas=thetas, endpoint = 'rate')
    power_summaries[[j]] = pp1; j = j+1;
    power_summaries[[j]] = pp2; j = j+1;
  }
  save(power_summaries, file = out_path)
} else {
  load(out_path)
}

Fig 2

# compute average 1 - type 2 error
power_list = lapply(power_summaries,
                    function(ll) apply(ll, 1:2, function(x) mean(x<my_alpha)))

cols=c(RColorBrewer::brewer.pal(8,'Dark2')[c(2,8)])
my_cols = rep(cols, length(Max_follows))
my_ltys = foreach(x = 1:length(Max_follows), .combine = c) %do% rep(x,2)


par(mfrow=c(2,3),las=1,lwd=1.5, cex.lab=1.6,cex.axis=1.4,
    family='serif',mar=c(6,5,3,3),bty='n')
for(i in 1:length(Effect_sizes)) {
  plot(Sample_sizes, power_list[[1]][i,], ylim = c(0,1),type='l',lwd=2,
       xlab='Sample size per arm',ylab='Power',col=my_cols[1], lty=my_ltys[1])
  for(j in 2:length(power_list)){
    lines(Sample_sizes, power_list[[j]][i,], col=my_cols[j], lwd=2,lty=my_ltys[j])
  }
  mtext(text = paste(100*(Effect_sizes[i]-1),'%',sep=''), side = 3, cex=1.5)
  if(i==1){
    legend('topleft',col=c(cols, rep('black',length(Max_follows))),bty='n',
           lty=c(1,1,1:length(Max_follows)),inset=0.03,title = 'Endpoint',
           legend = c('Time to clearance','Rate of clearance'),cex=1.3)
    legend(x = 20, y = 0.7,col='black',title = 'Duration of follow-up (days)',
           lty=1:length(Max_follows),inset=0.03,bty='n',
           legend = Max_follows, cex=1.3)
  }
}

Simulation 2 - how often to measure qPCR?

From twice daily to only twice over ten days for the estimation of the rate endpoint

# Fix effect size
Effect_sizes = 1.075
follow_days = lapply(c(0.5,1,2,5,10), function(x) seq(0,10,by=x))
power_rates = list()
out_path = 'Rout/frequency_follow_up.RData'
if(FORCE_RE_RUN | !file.exists(out_path)){
  for(i in 1:length(follow_days)){
    writeLines(sprintf('Doing simulation with %s follow-up points',
                       length(follow_days[[i]])))
    power_rates[[i]] = compute_pvalues(Effect_sizes = Effect_sizes,
                                       Sample_sizes = Sample_sizes,
                                       Follow_up_days = follow_days[[i]],
                                       Nsim = Nsim,
                                       summary_function = compute_rate,
                                       thetas=thetas, endpoint = 'rate')[1,,]
  }
  save(power_rates, file = out_path)
} else {
  load(out_path)
}

Plot output

summaries = lapply(power_rates, 
                   function(x) apply(x,1,function(x) mean(x<my_alpha)))

par(mfrow=c(1,1),las=1,cex.lab=1.5,cex.axis=1.5,
    family='serif',mar=c(5,5,3,2),bty='n')
mycols = RColorBrewer::brewer.pal(length(summaries), name = 'Dark2')
plot(Sample_sizes, summaries[[1]],type='l', col=mycols[1],panel.first=grid(),
     ylab='Power', xlab= 'Sample size per arm', ylim=c(0,1),lwd=3)
for(i in 2:length(summaries)){
  lines(Sample_sizes, summaries[[i]],col=mycols[i],lwd=3)
}
legend('topleft',legend = sapply(follow_days,length),col=mycols,cex=1.5,bty='n',
       inset=0.03,lwd=3, title='Number of viral load measurements\n taken over 10 days')

Simulation 3 - duration of follow-up

Fix effect and sample size and look at how duration influences power

Effect_sizes = 1.075
Sample_sizes = 50
Max_follow = seq(3, 21, by=1)
out_path = 'Rout/duration.RData'
Nsim = 10000
if(FORCE_RE_RUN | !file.exists(out_path)){
  power_TE_follow = power_rates_follow = array(dim = c(length(Max_follow), Nsim))
  # do for rate endpoint
  for(mm in 1:length(Max_follow)){
    power_rates_follow[mm,]=
      compute_pvalues(Effect_sizes = Effect_sizes,
                      Sample_sizes = Sample_sizes,
                      Follow_up_days = 0:Max_follow[mm],
                      Nsim = Nsim,
                      summary_function = compute_rate,
                      thetas=thetas, endpoint = 'rate')[1,1,,drop=T]
  }
  
  # do for time to event endpoint
  for(mm in 5:length(Max_follow)){
    power_TE_follow[mm,]=
      compute_pvalues(Effect_sizes = Effect_sizes,
                      Sample_sizes = Sample_sizes,
                      Follow_up_days = 0:Max_follow[mm],
                      Nsim = Nsim,
                      summary_function = compute_clearance_time,
                      thetas=thetas, endpoint = 'timetoevent')[1,1,,drop=T]
  }
  
  save(power_rates_follow, power_TE_follow, file = out_path)
} else {
  load(file = out_path)
}

Fig 3

par(mfrow=c(1,1),las=1,lwd=1.5, cex.lab=1.5,cex.axis=1.5,
    family='serif',mar=c(5,5,3,2),bty='n')
plot(Max_follow, apply(power_rates_follow, 1, function(x) mean(x<my_alpha)),
     type='l', ylim = c(0,1), ylab='Power',panel.first=grid(),
     xlab= 'Duration of follow-up (days)',col=cols[2], lwd=3)
lines(Max_follow,apply(power_TE_follow, 1, function(x) mean(x<my_alpha)),
      lwd = 3, col=cols[1])
abline(h=c(0.6,.44),v=c(10,12),lty=2)

legend('topright',legend = c('Time-to-clearance', 'Rate-of-clearance'),
       col=cols,lwd=3, title = 'Endpoint', cex=1.5)